Hybrid photon–phonon blockade

We describe a novel type of blockade in a hybrid mode generated by linear coupling of photonic and phononic modes. We refer to this effect as hybrid photon–phonon blockade and show how it can be generated and detected in a driven nonlinear optomechanical superconducting system. Thus, we study boson-number correlations in the photon, phonon, and hybrid modes in linearly coupled microwave and mechanical resonators with a superconducting qubit inserted in one of them. We find such system parameters for which we observe eight types of different combinations of either blockade or tunnelling effects (defined via the sub- and super-Poissonian statistics, respectively) for photons, phonons, and hybrid bosons. In particular, we find that the hybrid photon–phonon blockade can be generated by mixing the photonic and phononic modes which do not exhibit blockade.

www.nature.com/scientificreports/ mechanical systems 47 and studying single-phonon optomechanics, in addition to offering a source of single-or multiple phonons 50,51 . PB can be changed into light transmission 52 , e.g., by photon-induced tunnelling (PIT) 6 . This is another nonclassical photon-number correlation phenomenon, in which the probability of observing more photons in a higher manifold of the system increases with the generation of the first photon near the resonance frequency of the system. Multi-PIT effects were also predicted 29 , including those generated by squeezing 53 .
For simplicity, we use here the abbreviation PB, when referring to the blockade of not only photons, but also of phonons or hybrid photon-phonon bosons. The precise meaning can be found from its context, e.g., when we refer to a specific mode, including the optical (a), mechanical (b), or hybrid (c) modes. Analogously, PIT denotes a given particle-induced tunnelling among the three types of excitations.
Nanomechanical resonators can coherently interact with electromagnetic radiation 54 , and quantum correlations between single photons and single phonons were studied for a single entangled photon-phonon pair 55 or via photon and phonon blockade effects in optomechanical systems 56 . A mechanical switch between PB and PIT has been studied recently 57 . PB and PIT effects in systems comprising mechanical and optical resonators, which are characterised by the same or similar bare frequencies, to our knowledge, have not been studied experimentally yet, although they seem to be experimentally feasible and, thus, they are at focus of this paper.
Crucial signatures of PB and PIT can be observed by measuring the second-order correlation function, g (2) (0) . Specifically for photons, (1) the condition of g (2) (0) < 1 defines the sub-Poissonian photon-number statistics (also referred to as zero-delay-time photon antibunching), which indicates the possibility of observing PB, while (2) the condition g (2) (0) > 1 , defines the super-Poissonian statistics (also referred to as zero-delaytime photon bunching), which is a signature of PIT in a given system. To observe the 'true' effects of PB and PIT, also other criteria should be satisfied, such as nonzero-delay-time photon antibunching and higher-order sub-Poissonian photon-number statistics. Indeed, an ideal conventional PB, which can be served as a singlephoton source, usually should also be verified by studying higher-order correlation functions, g (n) (0) for n > 2 . For example, in case of single-PB (1PB) conditions g (2) (0) < 1 and g (n) (0) < 1 for n > 2 should be fulfilled.
PB can be verified also in other ways via demonstrating, e.g., a staircase-like dependence of the mean photon number (or measured power transmitted through a nonlinear resonator) on the energy spectrum of the photons incident on the resonator 8,52 . Such a dependence is the photon analogue of the Coulomb staircase. All of the above criteria are just necessary but not sufficient conditions for demonstrating PB. A sufficient condition could be, e.g., showing a high fidelity of a given generated light (with a nonzero mean photon number) to an ideally truncated two-dimensional state, which is the closest to the generated one. This approach was applied in, e.g., 26,35,36 . The latter two types of PB tests are, however, are not applied in this paper.
Conventional single-PB prevents the absorption of a second photon with a specific frequency due to the nonlinearity of a given system. Such a nonlinearity can be described by a Kerr-type interaction and/or can induced by an atom (real or artificial) coupled to a resonator. An artificial atom can be realised by, e.g., a quantum dot 23,58,59 in cavity QED 10 or a superconducting qubit or qudit in circuit QED 52 .
Unconventional PB, which is induced by destructive interference, operates better for very low (or even extremely low) mean photon numbers 12,13 . This can be disadvantageous by considerably decreasing the probability of generating a single photon. But, at the same time, it can be an advantage, because a very small mean photon number usually reduces the chance of generating multi-photon states and inducing higher-order coherence. This is not always the case, and even if the probability of observing two photons is suppressed, higher-order coherence might be enhanced, leading to the generation of multi-photon states 46 .
In this paper, we consider an optomechanical system, which generates photonic and phononic modes. Then we apply a balanced linear coupling transformation to the these modes to create hybrid modes (also referred to as supermodes). We study the interplay between photons and phonons resulting in their nonclassical number correlation effects. Thus, we find such system parameters to observe either PB or PIT in the four modes. In particular, we predict PB in one of the hybrid modes, but not in the individual (photon and phonon) modes, i.e., this PB is created from the two modes, which do not exhibit PB. We refer to this effect as hybrid photon-phonon blockade, which is the main result reported here.
Specifically, we define hybrid photon-phonon blockade as the blockade of hybrid-mode bosons (polaritons) obtained by coupling photons of an optical or microwave mode with phonons of a mechanical mode by a balanced linear coupler. The idea and criteria for testing this type of blockade are analogous to those for other known blockade effects (e.g., of photons, phonons, or magnons), but it is predicted for another type of bosons. We show that this hybrid blockade can occur by coupling the modes, which exhibit neither photon blockade nor phonon blockade.
To show this effect we analyse the system of two linearly-coupled resonators: a superconducting microwave resonator (SMR), which might be a transmission line resonator, and a micromechanical resonator, referred to as a quantum drum (QD), which is capacitively coupled to the SMR. To generate any kind of PB (including unconventional PB), one needs to incorporate a nonlinearity into a given system 17,60,61 . This can be done by coupling one of the resonators (e.g., the SMR) to a qubit (e.g., an artificial superconducting two-level atom). We also assume that the system is driven either at the QD or the SMR as described in detail in the next section.
The paper is organised as follows: first, the hybrid optomechanical system and its Hamiltonians are introduced. We also define the hybrid photon-phonon modes, which can be generated by the balanced linear coupling of photonic and phononic modes. Then, we study the correlation effects in the photonic, phononic, and one of the hybrid modes in the system driven at either the optical or mechanical resonator, respectively, for experimentally feasible parameters specified in "Methods". We then predict and analytically explain the generation of unconventional hybrid-mode blockade via a non-Hermitian Hamiltonian method. We systematically study different weaker and stronger criteria for observing blockade and tunnelling effects in our system. We also find all the eight combinations of the conventional blockade and tunnelling effects in the three modes. In particular, Figure 1 shows the schematics of the studied hybrid system, which consists of a superconducting two-level artificial atom (a qubit) embedded in a waveguide and coupled to an SMR, which might be a transmission-line resonator. This qubit induces anharmonicity in the SMR, which is crucial for observing PB. Our setup includes also a microwave-frequency mechanical resonator (a QD), which is capacitively coupled to the SMR. The nonlinearity of the QD is induced indirectly by the linear coupling of the QD to the effectively nonlinear SMR.

The system and Hamiltonians
The free Hamiltonian of the SMR is H a = ω SMR a † a , where ω SMR is its resonance frequency (assumed here of the order of tens of GHz) and a (a † ) is the photon annihilation (creation) operator. We can reasonably assume the SMR quality factor as Q SMR ≈ 10 4 . The free Hamiltonian of the QD is H b = ω m b † b , where ω m is its resonance frequency and b (b † ) is the phonon annihilation (creation) operator. In our numerical simulations, we set ω m /2π = 7.8 GHz and the QD quality factor as Q m ≈ 260 . Moreover a two-level quantum system has the ground state |g� and the excited state |e� with transition frequency ω q (set here of the order of ω m and ω SMR ). The free qubit Hamiltonian is described as H q = ω q σ + σ − , where σ + = |e��g| ( σ − = |g��e| ) is the atomic raising (lowering) operator. Thus, the total free Hamiltonian of the system is H 0 = H a + H b + H q . The complete Hamiltonian (without driving) of our coupled system can be given by ( = 1) which includes the three coupling terms: (1) the Jaynes-Cummings term describing the interaction between the SMR and qubit under the rotating-wave approximation (RWA); (2) the radiation-pressure term with coupling strength g r ; and (3) the Hopfield-type nonlinear coupling term with strength g l . The g l coupling can be realized via a capacitor, as shown in Fig. 1 and explained in a more detail in Ref. 48 for a similar system. Note that the g l term describes canonical position-position (momentum-position) interactions, where g l is real (imaginary) for H ′ + ( H ′ − ). These interactions can be interchanged by adding the π/2 phase to a, a † , i.e., a → ia and a † → −ia † . This extra phase does not change number correlations in the modes a and b. In typical ranges of parameters of analogous superconducting circuits 62 , the g l term is dominant, so the radiation-pressure term can be neglected 63 . Moreover, although the counter-rotating terms ab ± b † a † , which appear in the g l -interaction, play an important role in the ultrastrong and deep-strong coupling regimes 64 , but they can be safely omitted under the RWA, which is valid in the weak and strong coupling regimes. Indeed, the latter regimes are solely studied in this paper, as discussed below. Then the Hopfield nonlinear g l -interaction becomes effectively linearised. Thus, Hamiltonian (1) reduces to Figure 1. Schematics of the discussed circuit-QED-based realisation of the considered hybrid optomechanical system. It consists of a superconducting qubit embedded in a superconducting microwave resonator (SMR), e.g., a transmission-line resonator, to induce its nonlinearity. A quantum micromechanical resonator, which is referred to as a quantum drum (QD), is coupled to the SMR with a tunable capacitor C g . We assume that the system is driven either at the SMR or QD. Dashed semicircular curves visualise that the QD is oscillating. The driving and motion detection of the QD can be realised by controlling the static magnetic field B, potential V g , and alternating current I(t), as described in Ref. 47 for detecting phonon blockade. www.nature.com/scientificreports/ where the linear-coupling strength is denoted by f, which replaces the symbol g l . Analogously to g l , f is real (imaginary) for H + ( H − ). In the following, for simplicity, we focus on studying the canonical position-position interactions between the modes a and b, as described by H + . The eigenstates of Hamiltonian H ± can be referred to as atomic-optomechanical polaritons or atom-cavity-mechanics polaritons 65 . It is clear that Hamiltonians H ± conserve the polariton number, which is the total number of excitations. Thus, H ± can be diagonalised in each subspace (or manifold) H (n) with exactly n polaritons. The RWA is fully justified assuming both (1) the weak-or strong-couplings and (2) small detunings between the SMR and QD, and the SMR and qubit (see, e.g., Ref. 66 ). We stress that these conditions are fully satisfied for the parameters applied in all our numerical calculations in this paper. Thus, the Jaynes-Cummings and frequency-converter (or linear-coupler) models can be applied. However, the RWA cannot be applied in the ultrastrong and deep-strong coupling regimes, as defined by g > 0.1 ω i and g > ω i , respectively 64 , where i = SMR, m, q . In these regimes, the quantum Rabi and Hopfield models cannot be reduced to the Jaynes-Cummings and frequency-converter models, respectively. However, we study the system for the parameters specified in Eqs. (28)- (30), for which the ratios of the coupling strengths and frequencies, f /ω i and g/ω i , are < 0.002 . So, the system is in the strong-coupling regime, and far away from the border line with the USC regime. Moreover, the chosen detunings are |ω SMR − ω m |/ω SMR ≤ 2.6 × 10 −3 and |ω SMR − ω q |/ω SMR < 8 × 10 −4 . Thus, it is clearly seen that we can safely apply the RWA. Anyway, as a double test, we have calculated time-dependent secondorder correlation functions for the Hamiltonian H ′ ± and H ± for the parameters set in Eqs. (28)-(30) for various evolution times assuming classical drives (as specified below) and no dissipation. And we have found that the differences between the correlation functions calculated for the models with and without the RWA are negligible on the scale of figures. The inclusion of dissipation in the system makes such differences even smaller.
We assume that an optical pump field of frequency ω p is applied either to the SMR mode a, as described by or to the QD mode b, as given by to drive (excite) the system (with coupling strength η a or η b ) from its ground state and to induce the emission of photons and phonons. Thus, the total Hamiltonian becomes Direct driving of the QD can be implemented by a weak-oscillating current, as considered in Refs. 47,48 , where the drive strength η b is proportional to the current amplitude I(t) and the magnetic field B shown in Fig. 1. The SMR can be driven in circuit-QED systems in various ways 62 . Note that by driving directly the SMR (or alternatively the QD), one also indirectly drives the QD (SMR) through the capacitive coupling C g , as shown in the scheme in Fig. 1. So, by referring to the SMR-or QD-driven systems, we indicate only the resonator, which is directly pumped, although finally both resonators are driven.
The inclusion of an additional nonlinearity in the QD and/or applying drives to the qubit(s) and both resonators is not essential for the prediction of hybrid blockade, but this could enable achieving stronger photon-phonon antibunching and more sub-Poissonian statistics.
Considering the case, where the pump field drives only the SMR, to remove the time dependence of the Hamiltonian H (n) (t) and to obtain its steady-state solution, we transform the system Hamiltonian into a reference frame rotating at frequency ω p .
We apply the unitary transformation U R (t) = exp −iN polariton ω p t to H (n) according to the general formula Thus, H (a) (t) reduces the time-independent SMR-driven Hamiltonian: where � i = ω i − ω p for i = a, b, q . So, in particular, b ≡ m ( a ≡ SMR ) is the mechanical (microwave) resonator frequency detuning with respect to the pump frequency. Analogously, in the same rotating frame, H (b) (t) reduces to the QD-driven Hamiltonian:   www.nature.com/scientificreports/ We recall that Eqs. (8) and (9) are directly derived from Eq. (6) for H + given in Eq. (2). Moreover, H + is derived from Eq. (1) assuming the RWA, which is justified for small detunings in the weak-and strong-coupling regimes, which are the only numerically studied regimes in this paper, as emphasised above. Indeed, the studied ranges of parameters guarantee the system evolution is far from the USC regime. Note that the Hamiltonian H ′ + in (1) for g r = 0 with an additional drive term H rot given in Eqs. (8) and (9) but with the additional term f [ab exp(2iω p t) + h.c.] . In all our numerical calculations we set ω p of the order of GHz. Thus, the effect of this rapidly oscillating term is negligible compared to all the other terms in the Hamiltonians. Moreover, we have also assumed that the optomechanical term g r is negligible. In general, this assumption is not necessary, because the g r term can be reduced (in the red-detuned regime) to an interaction term describing a linear coupler (or a beam splitter), which can be combined with the f term. Anyway, for simplicity concerning both theory and potential experiments, we set g r = 0 . We have also assumed that the system is driven at either the mechanical or optical mode to obtain effectively time-independent Hamiltonians in a rotating frame. This simplification would not be directly possible by considering the system driven simultaneously at both modes with different frequencies. Figure 2 shows the structure of the energy spectrum for the hybrid system Hamiltonian (2). To study the sub-Poissonian light generation in hybrid modes, we apply to the SMR and QD modes a balanced linear coupling transformation, which is formally equivalent to a balanced (50/50) beam splitter (BS). This transformation creates the hybrid (or cross) photon-phonon modes:  www.nature.com/scientificreports/ The dynamics of an open system in the presence of losses under the Markov approximation can be described within the Lindblad approach for a system reduced density matrix ρ satisfying the standard master equation, which is given in terms of the Lindblad superoperator D where κ a , κ b , and γ are the decay rates for the SMR, QD, and qubit, respectively.
All our numerical calculations and their analyses are given for the system parameters, which satisfy the conditions for the weak or strong-coupling regimes and for small detunings between the SMR, QD, and qubit. Thus, we can safely apply the standard master equation given in Eq. (12). Of course, if one considers Eq. (1) for the system in the USC or deep-coupling regimes, then the master equation in Eq. (12), should be replaced by a generalised one, e.g., of Refs. 64,[67][68][69] .
We also note that the application even of a single classical drive to the Jaynes-Cummings model in the strong-coupling regime effectively creates counter-rotating terms, which can induce a variety of USC effects, as shown explicitly in Ref. 70 . Thus, to confirm the validity of our results, we have applied the generalised formalism described in Ref. 67 , which is valid for arbitrary light-matter coupling regimes, including the weak-, strong-, and USC regimes. In particular, we calculated the correlation functions g (n) (0) defined in terms of the positive-( X + n ) and negative-[ X − n = (X + n ) † ] frequency components of the canonical position operators: X a = a + a † for photons, X b = b + b † for phonons, and X c = c + c † for hybrid-mode bosons in the qubit-SMR-QD dressed basis. We calculated the steady states of the system by solving numerically the generalised master equation of Ref. 67 for the Hamiltonians H ′ and H ′′ . As expected from general considerations, our numerical calculations for the parameters set in Eqs. (28)-(30) using the standard and generalised formalisms based on H ′ (as well as H ′′ ) give effectively the same results.
In our simulations, we assume that the system is prepared in the ground state |n = 0, g�|m = 0� (i.e., with no photons in the SMR, no phonons in the QD, and the qubit is in the ground state), such that a given pump laser can drive the SMR photons in the microwave frequency range. Note that the choice of initial states affects the short-time evolution of our system, but has no effect on the steady-state solutions in the time limit, assuming the single-photon and single-phonon damping channels, as described in Eq. (12). However, as shown in Ref. 32 , initial states of a system can indeed affect steady states of the system, thus can also change PB, in case of quantum engineered dissipation channels allowing for, e.g., two-photon dissipation only.
In the following sections, we show that it is possible to observe both PB and PIT in the hybrid mode in the weak, mediate, and strong coupling regimes compared to the decay rates of the SMR, QD, and qubit. In particular, we show that the system can generate the hybrid photon-phonon modes with strongly sub-Poissonian (or super-Poissonian) statistics by mixing the SMR and QD modes with strongly super-Poissonian (or sub-Poissonian) statistics.

Hybrid-mode blockade in the SMR-driven system
Here we analyse in detail various blockade and PIT effects in the SMR-driven dissipative system described by the Hamiltonian H ′ and the master equation (12) for the parameters specified in Eq. (28).
Photon/phonon-number statistics of the modes generated by our hybrid system can be described quantitatively by calculating the zero-delay-time kth-order correlation function ( kth-order intensity autocorrelation function), where z = a, b, c, d and k = 2, 3, . . . . In the special case of k = 2 , which is of particular interest in testing single-PB and single-PIT, the three different types of the boson-number statistics can be considered: the Poissonian [if g (2) , and sub-Poissonian (otherwise). Analogously, one can define higher-order Poissonian, sub-Poissonian, and super-Poissonian statistics for k > 2 . Such higher-order criteria are not only crucial in analysing multi-PB and multi-PIT effects 11,29,53 , but they are also important in testing whether a specific PB effect is a 'true' PB, which can be used for generating single photons or phonons. These higher-order statistics are studied in "Methods". Figure 3(a) shows g (2) (0) as a function of the qubit-SMR coupling for the SMR-driven system with the parameters specified in Eq. (28). The regions, when the sub-Poissonian statistics in the hybrid mode c is accompanied by the super-Poissonian statistics in the modes a and b, are indicated by the yellow background in this and other figures. This area in yellow colour is referred to as Case 7 in Table 1, in which we observe strongly super-Poissonian photons (phonons) in the SMR (QD); whereas a single excitation is observed in the hybrid mode. The system parameters, which lead to Case 7, are found by numerical simulations and are discussed below.
Note that Fig. 3a shows these effects in the strong coupling regime 64 , i.e., when the qubit-SMR coupling constant g is larger than the system damping rates: g/κ max > 1 , where κ max = max{κ a , κ b , γ } . On the other hand, Fig. 3b shows the same yellow region in the weak-coupling regime, i.e., when g/κ max < 1 , but this figure was calculated for the QD-driven system, which is discussed in the next section.
By considering the values of Eq. (28), the SMR decay rate is κ a = 1.5γ , given that the mode a is always in the strong qubit-SMR coupling regime in the region of our interest. This results in Rabi-type oscillations of g (2) (0) that occur in the SMR mode a and the hybrid mode c. In Fig. 3a, both weak and strong coupling regimes are shown corresponding to g smaller or larger than the maximum decay rate of the whole system. www.nature.com/scientificreports/ Given the set of parameters in Eq. (28), we are in the good-cavity regime 71 , because κ a < {κ b , g, f } . In the range g/2π ∈ (4.5, 42) MHz, the hybrid mode c has the sub-Poissonian statistics, while the SMR mode has the super-Poissonian statistics in all the shown cases and a very weak sub-Poissonian statistics occur for phonons in the QD mode b, but still corresponding to Case 4 in Table 1. This behaviour changes to the super-Poissonian statistics in the mode b, which corresponds to Case 7, as shown in Fig. 3a. There is a transition for the mode c from the sub-Poissonian to super-Poissonian statistics, which corresponds to switching from Case 7 to Case 8 in the strong-coupling regime, where the other two modes are both super-Poissonian. Observing g (2) (0) > 1 witnesses PIT and the quantum nature of this effect is explored further below.
In order to better probe and understand the dynamics of the system in specific parameter regimes, we analyse also the delay-time second-order photon correlation function, defined as where n z (t) = z † (t)z(t) is the boson number in the modes z = a, b, c, d , and the operator products are written in normal order (::) and in time order T . With g (2) z (τ ) another quantum optical number-correlation phenomenon can be investigated. Specifically, in case of photons, it is referred to as photon antibunching if g (2) (0) < g (2) (τ ) , photon unbunching if g (2) (0) ≈ g (2) (τ ) , and photon bunching if g (2) (0) > g (2) (τ ) , which is usually defined for short or very short delay times τ 72 . It is worth noting that photon antibunching was first experimentally observed in the 1970s by Kimble, Dagenais, and Mandel 73 . This was historically the first experimental demonstration of the quantum nature of an electromagnetic field, which cannot be explained classically, unlike photoelectric bunching.
Analogously, one can also investigate the antibunching and bunching of phonons and/or hybrid-mode bosons. Note that the term photon antibunching is often interchangeably used with the sub-Poissonian photon-number i (0) (in the common logarithmic scale) versus the ratio of qubit-SMR coupling strength and the largest decay rate. Different predictions of the sub-and super-Poissonian boson number statistics, which can be interpreted, respectively, as the PB and PIT effects, of the photonic (a), phononic (b), and hybrid (c) modes assuming: (a) the SMR-driven system with parameters specified in Eq. (28) and (b) the QD-driven system with Eq. (29). All the shown Cases (i.e., 4, 6, 7, and 8) correspond to those listed in Table 1. The broken line at g = max κ j is the border line between the strong-and weak-coupling regimes. Table 1. Different predictions of the super-and sub-Poissonian particle (i.e., photon, phonon or hybrid photon-phonon)-number statistics (PNS) corresponding, respectively, to PIT and PB, for the photon mode a, phonon mode b, and hybrid photon-phonon mode c, where c (0) − 1] and the last column indicates each prediction of the mode a, b, and c in the specific colour that is used in our plots. All these cases can be seen in Fig. 10. www.nature.com/scientificreports/ statistics 21 . However, to avoid confusion, one can refer to single-time (or zero-delay-time) photon antibunching if defined by g (2) (0) and two-time (or delay-time) photon antibunching if defined via g (2) (τ ).
In Fig. 4, we plotted g (2) (τ ) for the range [0, 1.5] of g/κ max . This range is also shown in Fig. 3a, where the examples of Cases 4 and 7 can be identified. As expected, one can see oscillations in the SMR and hybrid modes in Figs. 4a,c, respectively. These oscillations are induced by the competition between the qubit-SMR coupling g and the SMR-QD hopping f in our system. Apparently, by analysing g (2) (τ ) in the weak-coupling regime, the frequency of the oscillations is smaller than that in the strong-coupling regime, in which the oscillations are caused by both couplings g and f. Moreover in a very weak coupling regime, where g ≪ 1 oscillations occur due to the hopping strength f, with the period 2π/f 74 . This means that, in the weak-coupling regime, also the coupling between the SMR and QD can generate oscillations in our system, where in this case the period of oscillations, which are induced by f = 5.5γ , is approximately equal to τ ≈ 0.036 , which coincides with the period deduced from the graph, as seen in Fig. 5c. These detrimental oscillations should be suppressed on a time scale longer than the SMR lifetime τ = 1/κ a to enable boson antibunching to survive in the area of our interest.
Various combinations of correlations effects are shown in Fig. 5. All panels in Fig. 5 show that the photon mode a is super-Poissonian and bunched, while the hybrid mode c is sub-Poissonian and antibunched. However, the properties of the phonon mode b are different in every panel. Specifically, the mode b is in panel: a (τ ) for the photonic mode, (b) g (2) b (τ ) for the phononic mode, and (c) g (2) c (τ ) for the hybrid mode versus the coupling strength g and the delay time τ . We consider here the SMR-driven system with parameters specified in Eq. (28), which enable us to observe the single-photon resonances in the mode c. For clarity, all the values of the correlation functions ≥ 2 are truncated at 2.  i (τ ) for the SMR mode a, the QD mode b, and the hybrid mode c modes assuming: (a-c) the SMR-driven system specified in Eq. (28) with f = 5.5γ and κ max = κ b = 6γ , and (d,e) the QD-driven system in Eq. (29) with κ max = 7.5γ , where we additionally set: b (τ ) for non-zero but short delay times τ ], (b) Poissonian and unbunched, (c) sub-Poissonian and unbunched, (d) super-Poissonian and bunched, and (e) Poissonian and bunched, as usually considered for very short delay times τ . Note that panels (a, b, c) are for the SMR-driven system, while the remaining panels (d, e) are for the QD-driven system, which are discussed in detail in the next section.
In particular, it is seen that by decreasing the coupling at g/κ b = 1.1 in Fig. 5b, the QD mode b is unbunched with the Poissonian statistics, while the hybrid mode c exhibits antibunching g (2) (0) < g (2) (τ ) and the sub-Poissonian statistics g (2) (0) < 1 , in both cases. The role of the auxiliary mode b is, in a sense, to convert the super-Poissonian into sub-Poissonian statistics in the mode c.
The destructive interference of both modes a and b, at the balanced linear coupler, can result in the sub-Poissonian statistics of the hybrid modes. We observe this effect even in the weak-nonlinearity (or weak-coupling) regime, which witnesses unconventional PB, as discussed in detail in "Methods". It is worth noting that in this study we are aiming at observing g (2) (τ ) < 1 not only at τ = 0 , but also for non-zero delay times (e.g., τ ∈ [0, 0.1] ), as in standard experimental demonstrations of the boson antibunching statistics reported in, e.g., Refs 7,75 . Thus, the cases shown in Fig. 4a,c can hardly be considered as convincing demonstrations of the sub-Poissonian statistics, because of the oscillations, which occur in g (2) a,c (τ ) with increasing τ . More convincing demonstrations of these effects without such oscillations (or by considerably suppressing them) are presented in Figs. 6 and 7, as analysed in detail in the next section.
To explain the super-Poissonian photon-number statistics and photon bunching in the mode a for the system pumped in the SMR mode, let us analyse Fig. 5a with g ≈ κ m concerning the anharmonicity of the energy levels in these cases.
The g term in Eq. (2) corresponds to the standard Jaynes-Cummings model with the familiar eigenvalues 62 : with the corresponding eigenstates: which are often referred to as dressed states or dressed-state dublets, where θ n = � n /� 1 is the mixing angle, � 1 = ω q − ω SMR is the detuning between the SMR and qubit. Moreover, n = 2g √ n + 1 can be interpreted as the n-photon Rabi frequency on resonance, so, in particular, 0 = 2g is the vacuum Rabi frequency. Thus, the energy spectrum is clearly anharmonic, which is a necessary condition to observe PB. Note that the Jaynes-Cummings interaction can be effectively described in the dispersive limit (i.e., far off resonance) as a Kerr nonlinearity (for a detailed derivation see, e.g., 50 ), which is the standard nonlinearity assumed in many predictions of PB effects.
To demonstrate the anharmonic energy levels of the complete Hamiltonian H + on resonance (see Fig. 2), we assume a weak drive coupling strength η a . Given that, the system Hilbert space can be truncated. We assume that the polariton number is at most equal to two in this weak-drive regime. The ground state is |ψ 0 � = |0, 0, g� with the corresponding eigenvalue E 0 = 0 . The three eigenvalues of the first manifold (with eigenstates containing a single polariton), as shown in Fig. 2b, are: while the five eigenvalues of the second manifold (with eigenstates containing two polaritons), which are shown in Fig. 2c, read: 16) |n, +� ≡ cos θ n 2 |n�|e� + sin θ n 2 |n + 1�|g�, |n, −� ≡ − sin θ n 2 |n�|e� + cos θ n 2 |n + 1�|g�, Figure 6. Same as in Fig. 4, but for the QD-driven system with parameters given in Eq. (29). We observe here single-PRs and the corresponding single-PB effects. www.nature.com/scientificreports/ where f 1 = 3f 2 (10g 2 + 3f 2 ) + g 4 . In particular, by assuming f = 5γ and g = 7.5γ , the eigenenergies of the first and second manifolds are, respectively: (1) , � ± 9.01388γ ≈ � ± 9γ , and (2) 2 , 2� ± 5.82965γ ≈ 2� ± 6γ , and 2� ± 16.11725γ ≈ 2� ± 16γ. A simple way to probe the pumped mode is to record the second-order correlation g (2) (0) as a function of SMR , where the pump frequency ω p is changing (see Fig. 8). To do so, we first consider the resonance case as ω SMR = ω m = ω q = ω in Eq. (8) and ω − ω p = � . As depicted in Fig. 8a, one can see local minima with negative values in log g (2) (0) for the three modes, which indicate Case 1 in Table 1, at � SMR /γ = ±9 , which correspond to � = ± g 2 + f 2 ≈ ±9γ , given Eq. (17). This means that the pump frequency is located at the two dressed state dublets with energies E 1 and E 3 . And we are off-resonance from the second energy manifold, which implies the possibility of observing PB at these frequencies.
Furthermore, our simulations predict a maximum of log g (2) (0) ≈ 3 showing a strong super-Poissonian statistics in the three modes (corresponding to Case 8 in Table 1) as SMR → 0 . In particular, at � SMR /γ ≈ ±6 , the pump frequency is near E 4 ≈ −6 , respectively, of the second manifold, in which the probability of the two-photon resonance is maximised, as a signature of PIT. It signifies that the pump is in resonance with one i (τ ) (in the logarithmic scale) for the SMR mode a, the QD mode b, and the hybrid mode c modes in the QD-driven system assuming that g/κ max is equal to: (a) 3, (b) 2.1, (c) 3.8, and (d) 2.2. The four different predictions of correlations for the QD mode b correspond to all the cases listed in Table 2. (e,f) Same as in Fig. 4, but for the parameters given in Eq. (30). Note that panels (a-d) show the cross-sections of the 3D plot in (f) at the values of g/κ max marked by broken lines. www.nature.com/scientificreports/ of the levels in the second manifold of the hybrid system energy levels, here specifically E 1 and E 4 . One can see in Fig. 8, peaks (global maxima in the analysed range) of log g (2) n (0) > 0 for n = a, b, c at SMR = 0 . In particular, the probability of absorbing a single photon decreases here. However, if a photon is absorbed, it enhances the probability of capturing subsequent photons, this effect produces the super-Poissonian statistics, which is due to the fact that the probability of observing a single photon is also very small ( P 10g ≪ 1 ) and smaller than the probability of observing two photons 6,76 .
It is seen that, by tuning the drive frequency to the transition E 2 − E 0 in the energy spectrum of the total nonlinear system, the probability of admitting two photons increases. This results in the super-Poissonian statistics, which is opposite to the case, when the drive frequency is tuned to the transition E 1 − E 0 , when the probability of admitting subsequent photons decreases resulting in PB.
By assuming the off-resonance condition, ω SMR = ω m = ω q , we show in Fig. 8b the correlation functions for the three modes (a, b, c) as a function of SMR in the case, when the drive is tuned in-between the dressed state eigenenergies of the hybrid system.
The PB and PIT effects observed in Fig. 8 can be explained by considering some measures of the distances from resonances, as shown in Fig. 9a. The distances of the single-, two-, and three-photon resonances (PRs) are defined here, respectively, as: where ω p is the frequency of the pump that is tuned with respect to the energy of the hybrid system. Here ω (n) i are the frequencies (labelled with subscript i) in the nth manifold, so the minimalization is performed over ω (n) i for a given manifold n. Figure 9 shows the resonance distances versus SMR , where ω p is tuned with respect to the energy of the whole system. The dip in g (2) (0) at � SMR /γ = 10 (see Fig. 8b), which is characteristic for PB, corresponds to the resonance for a single excitation, as seen from D 1PR , and is off-resonance for higher excitations at that frequency (see Fig. 9a). The second-order correlation function g (2) c (0) for the hybrid mode has a dip as a signature of PB around � SMR /γ = −3.4 , while the modes a and b exhibit the super-Poissonian statistics (indicating PIT), as shown in Fig. 8b. This effect is witnessed as a dip in D 1PR and it is off-resonance for D 2PR and i (0) versus the frequency detuning SMR (in units of the qubit decay rate γ ) between the drive and SMR for: (a) the resonance case ω SMR = ω m = ω q (so also SMR = m = q ) and (b) the nonresonance case ω SMR = ω m = ω q , where ω b /γ = 1560 MHz. Note that by changing the pump frequency, different detunings appear with respect to the modes a and b, and qubit. We set g = 7.5γ and other parameters are given in Eq. (28). The numbering of the coloured regions correspond to the cases listed in Table 1. www.nature.com/scientificreports/ D 3PR , as illustrated in Fig. 9a, while the modes a and b exhibit PIT. This type of unconventional PB is discussed further in sections below. By decreasing � SMR /γ from 0 to − 2, the correlation function g (2) a (0) for the SMR mode in Fig. 8a resembles a shoulder in shape. We observe PIT at this point or region, as expected from our findings in the resonancedistant diagram in Fig. 9a. Indeed, there is a dip in D 2PR for higher resonances at this point, which explains the occurrence of PIT.
Let us consider now � SMR /γ → 3 in Fig. 8b for the pump frequency in resonance with the qubit, q = 0 , which is close to the resonance frequency of the hybrid mode. In this case multi-photon transitions are induced, which result in PIT at � SMR /γ = 3 , and we observe a peak in log g (2) (0) > 0 at this frequency in Fig. 8b. Clearly, we are here in resonance with higher-energy levels, while the drive strength is very small, η a /γ = 0.7 . The probability of observing a single photon is also small as the peak for c = 0 , but if a single photon is absorbed, then the probability of capturing subsequent photons increases, as for PIT.
The analysed system parameters are found by optimising our system to observe the super-Poissonian statistics in the SMR and QD modes. At the sub-Poissonian statistics area of g (2) (0) , it is possible to observe in Fig. 14 (in "Methods") that g (3) (0) > 1 and/or g (4) (0) > 1 , which are signatures of higher-order photon/phonon resonances and multi-PIT (see "Methods"). Actually, by calculating the second-order correlation function to witness the PB and PIT phenomena, higher-order correlation functions can be used to test whether a given effect is indeed: (1) single-PB or single-PIT, (2) multi-PB or multi-PIT, or (3) nonstandard versions of these effects, as discussed in "Methods" and, e.g., in Refs. 29,53 . As mentioned above, these parameters allow us to achieve the sub-Poissonian statistics for a relatively long delay times.

Hybrid-mode blockade in the QD-driven system
In this section, we analyse steady-state boson-correlation effects, including the hybrid-mode blockade and PIT, in the QD-driven dissipative system, as described by the Hamiltonian H ′′ and the master equation (12) for the parameters specified mostly in Eqs. (29) and (30).
To eliminate or at least to suppress the undesired oscillations in g (2) (τ ) , we assume in this section that our system is driven classically at the QD. Moreover, we assume that the SMR is in the bad-cavity regime, as κ SMR ≫ g 2 /κ SMR ≫ γ 71 . So, we apply the effective system Hamiltonian in the rotating frame, as given by Eq. (9). Even if the lifetime τ SMR = 1/κ SMR of the SMR is much shorter than that assumed in the SMR-driven system, which was discussed in the former section, the hybrid mode, as we show below, reveals no oscillations for quite long delay times, which is due to driving the QD.
To study boson-number statistics of our system, we compute the second-order correlation function g (2) (0) for the optimised parameters, which enables us to demonstrate Cases 4, 6, and 7 of Table 1 in Fig. 3b. In Case 7, which is of our special interest, the modes a and b are super-Poissonian, as log g (2) (0) > 0 , while the hybrid mode c is sub-Poissonian, as log g (2) c (0) < 0 . By increasing the coupling g between the SMR and qubit, the mode b becomes sub-Poissonian, as being affected by the nonlinearity of the mode a.
To check the second criterion for PB, the second-order correlation function g (2) (τ ) is considered below. Figure 6 shows g (2) (τ ) corresponding to g (2) (0) plotted in Fig. 3b showing Cases 4, 6, and 7. As expected, boson antibunching is observed for the hybrid mode, as shown in Fig. 6c, while the SMR mode reveals bunching, as illustrated in Fig. 6a. Moreover both phonon antibunching and bunching, in addition to unbunching [i.e., g (2) b (0) ≈ g (2) b (τ ) for τ 0 ], have been observed in the studied region of the QD mode, as shown in Fig. 6b. It is clear from Fig. 6 that the antibunching of bosons in the three modes survives in some specific coupling regime (around g = 0.7κ m ) for a relatively long delay time τ > 1/κ and oscillations in g (2) c (τ ) are absent in the hybrid mode c. Moreover, boson bunching is observed, when g (2) a (τ ) drops rapidly for delay times greater than the cavity photon lifetime, as considered in Fig. 5d,e.
To understand the delay-time dependence of the hybrid mode c, we consider Eq. (9), when the SMR, QD, and qubit have the same resonance frequency, ω SMR = ω m = ω q = ω and g = 4.5γ . As illustrated in Fig. 10a, there are three dips (local minima) in g (2) b (0) < 0 for the mode b of the QD, where we assumed g < min{κ a , κ b } and f > g . For these parameters, only a weak nonlinearity is induced in the mode b. Thus, the anharmonicity of energy levels cannot explain the PB effect observed as a dip at these three dips (see Fig. 9b). Actually, these dips in log g (2) b (0) are due to single-photon resonant transitions, which correspond to unconventional PB, as explained by the non-Hermitian effective Hamiltonian method in the next section and in "Methods". Figure 10c shows log g (2) i (0) for the three modes as a function of SMR . In this case, we assume that the resonance frequencies of the SMR, QD, and qubit are not the same, and the detuning of each mode with respect to ω p is different. It is shown that, when � SMR /γ → 2 , multiphoton transitions (and so PIT or multi-PB) can be induced in the mode a, where the pump frequency is in the resonance with the qubit, ω p = ω q . This effect is seen in Fig. 14 (in "Methods") corresponding to a local maximum in higher-order moments g (3) i (0) and g (4) i (0) . Likewise the resonance case, unconventional PB in the modes b and c can be explained by the method applied in the next section.
In Fig. 11, we study how the second-order correlation functions reveal the PIT regime, which corresponds to Case 8 in Table 1 ] for small pump strengths η a,b . By increasing the driving power at least to some values, which can be identified in the figures for specific modes, we observe that the correlation functions g (2) (0) also decrease for all the modes (except the mentioned case). This property confirms the nonclassicality of the predicted PIT in the hybrid system according to an additional criterion of 'true' PIT of Ref. 14

Unconventional blockade explanation via non-Hermitian Hamiltonian approach
In this section, we apply the analytical mathematical formalism of Ref. 45 , based on an non-Hermitian Hamiltonian, to identify the quantum interference effect that is responsible for inducing unconventional PB, i.e., strongly sub-Poissonian statistics in the weak-coupling regime or the weak-nonlinearity regime. We stress that this is an approximate approach, where the effect of quantum jumps is ignored 77,78 . By considering the system studied in the former section under the weak-pump condition, we can truncate the Hilbert spaces for the modes a and b and the qubit at their two excitations in total. This allows us to consider the total-system Hilbert space of dimension 3 × 3 × 2 = 18 . Moreover, the weak-pump condition implies that C 00g ≫ C 10g , C 01g , C 00e ≫ C 11g , C 10e , C 01e , C 20g , C 02g . Thus, the steady-state of the coupled system can be expressed as where |n a , n b , g/e� is the Fock state with n a photons in the SMR, n b phonons in the QD, and the lower ( |g� ) or upper ( |e� ) state of the qubit. The effective non-Hermitian Hamiltonian of the system can be written as where H ′′ is given by Eq. (9). Analogously, one can consider the non-Hermitian Hamiltonian with H ′ , given by Eq. (8).
In the weak-pump regime, the mean number of photons and phonons in the SMR and QD can be approximated as �n a � ≈ |C 10g | 2 and �n b � ≈ |C 01g | 2 , respectively. As derived in detail in "Methods", the second-order correlation functions for generated photons and phonons, under the same weak-pump conditions, can be given by: where the superposition coefficients C n,m,g are given in Eqs. (39) and (41). Figure 10. Correlation functions log g (2) i (0) versus the frequency detuning SMR (in units of the qubit decay rate γ ) between the drive and SMR for the QD-driven system for: (a,b) the resonant case with ω SMR = ω m = ω q (so also SMR = m = q ), and (c) the nonresonant case with ω SMR = ω m = ω q . Parameters are set in: Eq. (29) with g = 4.5γ for (a,c), and Eq. (30) with g = 9.5γ for (b). Eight different predictions, which correspond to all the cases listed in Table 1, are marked for the sub-and super-Poissonian number statistics in the photonic (a), phononic (b), and hybrid photon-phonon (c) modes. www.nature.com/scientificreports/ The hybrid photon-phonon modes, which are defined in Eq. (10), are the output modes of the balanced linear coupler with the SMR and QD modes at its inputs. As shown in "Methods", we find, analogously to Eq. (22), the second-order correlation function for the hybrid mode c reads: where the superposition coefficients C ′ n,m,e/g are given in Eqs. (40) and (41), and the sixth formula in Eq. (34). This approach enables us to explain unconventional PB generated in the hybrid system, which is the result of a destructive quantum interference effect that assures, together with other conditions, that the probability amplitude of having two photons in the SMR and QD is negligible. This method can also be used to find some optimal parameters to observe PB in the system. Figure 12 presents a comparison of our predictions based on the precise numerical solutions of the master equation in Eq. (12), as shown by thin curves, with those calculated from Eqs. (22) and (23) using the non-Hermitian Hamiltonian approach, as shown by thick curves. The locations of the maxima and minima of the correlation functions are found similar according to both formalisms. However, these extremal values can differ more distinctly, especially for the two global minima in the sub-Poissonian statistics of the mode b and the super-Poissonian maximum of the mode a. The differences result from the effect of quantum jumps, which are properly included in the master-equation approach and totally ignored in the non-Hermitian Hamiltonian approach (Fig. 12).

Different types of blockade and tunnelling effects
The sub-Poissonian statistics of a bosonic field, as described by g (2) (0) ≪ 1 , is not a sufficient criterion for observing a 'true' PB, which can be a good single-photon or single-phonon source. In fact, other criteria, such boson antibunching, g (2) (0) < g (2) (τ ) , and the sub-Poissonian statistics of higher-order correlation functions, g (n) (0) ≪ 1 , should also be satisfied (see "Methods"). Anyway, most of the studies of PB, and especially those on unconventional PB, are limited to testing the second-order sub-Poissonian statistics described by g (2) (0) < 1.
As explicitly discussed in Refs. 21,72,79,80 , photon antibunching and sub-Poissonian statistics are different photon-number correlation effects. So, the four cases listed in Table 2, can be considered as different types of PB and PIT. We show that all these effects can be observed in the studied system. For brevity, Table 2 is limited to phononic effects. PB, as defined in Case I and often referred to as a 'true' PB, can be a good single-photon sources; but, as mentioned above, other higher-order criteria should also be satisfied.
To show these four different effects, we use the parameters set in Eq. (30), where κ b ≪ κ a at the κ b = 0.002 γ , which indicates that the quality factor is Q ≈ 200, and so η b /κ b ≈ 100 in the case of a strong pump driving the QD mode with η b = 0.22 γ . Apart from the previously mentioned phenomena, such as observing the super-Poissonian statistics and bunching in the SMR and QD modes, while a hybrid mode exhibiting the sub-Poissonian Figure 11. Second-order correlation functions log g (2) i (0) versus the drive strengths: (a,c) η a for the SMRdriven system and (b,d) η b for the QD-driven system. Parameters are given in: (a) Eq. (28) with g = 7.5γ and ω p = 1554γ , which implies � q = −3γ , � b = 6γ , a = 0 ; (b) Eq. (29) with ω p = 1568γ , which implies q = 0 , � b = −8γ , and � a = 2γ ; (c) Eq. (28) with g = 7.5γ and ω p = 1551γ , which implies q = 0 , � b = 9γ , and � a = 3γ ; and (d) Eq. (29) with ω p = 1570γ , which implies � q = −2γ , � b = −10γ , and a = 0. www.nature.com/scientificreports/ statistics and boson antibunching, we find the four types of PB/PIT in the mode b in different coupling regimes, as shown in Table 2, which includes the examples of specific experimentally feasible values of g/κ a . Case I corresponds to a stronger form of PB, which we refer to as a 'true' PB, when the nonclassical nature of bosons is revealed by both their antibunching and sub-Poissonian statistics. Case II corresponds to a stronger form of PIT, which can be called a 'true' PIT, when bosons exhibit both classical effects: the super-Poissonian statistics and bunching. In Case III, one can talk about a weaker form of PIT or, equivalently, another weaker type of PB, as such bosons are characterised by the classical super-Poissonian statistics and their nonclassical nature is revealed by antibunching. Case IV represents another weaker form of PB or, equivalently, of PIT, which is characterised by the nonclassical sub-Poissonian statistics of classically bunched bosons. These results imply that one cannot say in general that the antibunching of bosons leads to their sub-Poissonian statistics and vice versa 21,79 .
Our focus in this paper is on the generation of PB in the hybrid mode, while the other two modes exhibit PIT. Note that this a very special case of Table 1, which shows that eight combinations of boson number correlation phenomena in the modes a, b, and c can be generated in our system, as specified by the numbered coloured regions in various figures corresponding to the cases in Table 1. Thus, we found all the eight possible combinations of the PIT and PB effects in the hybrid system for the parameters specified in Eqs. (28), (29), and (30).

Detection of the hybrid-mode correlation functions
Here, we describe two detection schemes for measuring the intensity autocorrelation functions for the hybrid photon-phonon modes c and d, as shown in Fig. 13.  i (0) versus the frequency detuning SMR in units of γ for the QD-driven system for the resonant case with ω SMR = ω m = ω q = γ × 1560 MHz (so also a = b = q ). The thin curves in each mode are obtained using the master equation in Eq. (12) and the thick curves are obtained from the non-Hermitian Hamiltonian method using Eqs. (22) and (23). Parameters are set in Eq. (29) except g = 4.5γ and κ a = κ b = 6γ. Table 2. Different single-and two-time phonon-number correlation effects induced in the QD mode, which can be observed for different values of the qubit-SMR coupling strength g with respect to the SMR decay rate κ a , e.g., by setting the other parameters to be the same as in Eq. (30). Here, PNS stands specifically for the phonon-number statistics of the mode b. Note that we also found examples of Cases I, II, and IV for the modes a and c using the same system parameters as for the mode b. Sub-Poissonian PNS g Phonon antibunching g Phonon bunching g (2) b (τ ) < g  www.nature.com/scientificreports/ The measurements of g 2 (τ ) for the photonic mode a and the phononic mode b are quite standard and are usually based on the Hanbury-Brown and Twiss (HBT) optical interferometry and its generalised version for phonons 81 , respectively. However, the measurement M (as schematically shown in Fig. 13a) of g (2) (τ ) , or even g (2) (0) , for the hybrid photonic-phononic modes c and d is quite challenging if applied directly. Here we propose two detection methods, as shown in Fig. 13b,c, for indirect measuring of g (2) c,d (0). The first operation of the measurement unit M in both schemes is a linear-coupler transformation of the hybrid modes (c, d) into (a ′ , b ′ ) , which, assuming that the process is perfect, should be equal to the original purely photonic (a) and phononic (b) modes.
We consider a linear coupler (formally equivalent to a beam splitter) described by a unitary operation U LC (θ) , which transforms the input operators a and b into: for a real parameter θ , where T = cos 2 θ and R = 1 − T = sin 2 θ are the transmission and reflection coefficients of the linear coupler, respectively. The studied hybrid modes are the special cases of Eq. (24) for c ≡ c(θ = π/4) and d ≡ d(θ = π/4) . Clearly, the first transformation U LC (−π/4) in Fig. 13b,c, is the transformation inverse to that in Fig. 13a.
Detection method 1 based on measuring photons and phonons. The correlation functions g (2) c,d (0) in the hybrid photon-phonon modes can be measured indirectly, as indicated in Fig. 13b, by measuring the observables: where k, l, m, n = 0, 1, 2 , by using the relations: and analogous relations for the hybrid mode d. The measurement units M ′ a and M ′ b in this method, as shown in Fig. 13b, describe the measurements of photons and phonons, respectively. It is seen that, in this approach, to determine g (2) c,d (0) , one has to measure the following observables: f 01 , f 10 , f 11 , f 02 , f 20 , f 12 , f 21 , and f 22 . Almost each observable f kl should be measured simultaneously with a specific observable g mn , which can be realised by a coincidence and count logic (CCL) unit in Fig. 13b.
The measurements of all the required photonic observables f kl can be performed by using, e.g., the Shchukin-Vogel method, which is based on balanced homodyne correlation measurements 82 . According to that method, a photonic signal is superimposed on a balanced beam splitter with a local oscillator, which is in a coherent state |α = |α| exp(φ)� with a tunable phase φ . A desired mean value of the observable f kl can be obtained by linear combinations of the coincidence counts registered by specific detectors for different local-oscillator phases φ . This part of the method corresponds to a Fourier transform. The simplest nontrivial configuration, which enables the measurement of the observables f 10 , f 01 , f 20 , and f 02 , requires four detectors and three balanced BSs, where additional input ports are left empty, i.e., allowing only for the quantum vacuum noise. By replacing the four detectors with four balanced BSs with altogether eight detectors at their outputs, one can measure any observable f kl for k + l ≤ 4 . These include the desired observables f 21 , f 12 , and f 22 . Of course, the observable f 22 can be measured in a simpler way via the HBT interferometry. The measurement of phononic observable g mn can be performed analogously just by replacing the balanced BSs by balanced phonon-mode linear couplers and using phonon detectors as, e.g., in Ref. 81 . The measurement of two-mode moments f kl g mn is, at least conceptually, a simple generalisation of the single-mode methods relying on proper coincidences in photonic and phononic detectors. Note that a multimode optical version of the original single-mode method was described in Ref. 83 . Figure 13c shows another realisation of the measurement unit M, to determine g (2) c,d (0) , and even g (2) c,d (τ ) . This method is, arguably, simpler and more effective than detection method 1, because it is based on measuring only photons and using standard HBT interferometry. Our approach was inspired by Ref. 48 , where the measurement of single-mode phonon blockade was described via an optical method instead of a magnetomotive technique, which was described in Ref. 47 , where phonon blockade was first predicted.

Detection method 2 based on measuring only photons.
Our measurement setup realises the following three transformations: (1) converting the phononic mode b ′ into a photonic mode b ′′ , (2) mixing the optical modes a ′ and b ′′ on a balanced BS to generate the modes c ′ and d ′ , which, in an ideal case, have the same boson-number statistics as the original hybrid photon-phonon modes c and d; and finally, (3) applying the conventional optical HBT interferometry for these two optical modes. In unit (1), this conversion corresponds to a multi-level SWAP gate, which can be implemented by a photonic-phononic linear coupler for θ = π/2, assuming that the auxiliary input mode e ′ is in the photonic vacuum state, while the output mode e ′′ is in the phononic vacuum state. In unit (3), the balanced BS action on the optical modes a ′ and �c † c� = 1 2 �f 11 � + �g 11 � + �f 01 g 10 � + �f 10 g 01 � , �c †2 c 2 � = 1 4 �f 22 � + 4�f 11 g 11 � + �g 22 � + 2�f 01 g 21 � + 2�f 10 g 12 � + �f 20 g 02 � + �f 02 g 20 � + 2�f 21 g 01 � + 2�f 12 g 10 � , www.nature.com/scientificreports/ b ′′ in Fig. 13c corresponds to the transformation of the balanced linear coupler on the photonic (a) and phononic (b) modes, as shown in Fig. 13a. Clearly, the linear-coupler transformation U LC (θ) is applied not only to the modes (a, b), but also to other modes. Thus, Eq. (24) should be adequately modified by replacing (a, b) by (c, d), (e ′ , b ′ ) , and (a ′ , b ′′ ) . For brevity, we omit their explicit obvious definitions here. Note that U LC (π/2) and U LC (π/4) correspond to a multi-level SWAP and Hadamard-like gates, respectively; while the balanced BS in Fig. 13c corresponds to U LC (π/4).

Discussion
We proposed a novel type of boson blockade, as referred to as hybrid photon-phonon blockade, which is a generalisation of the standard photon and phonon blockade effects. We predicted the new effect in a hybrid mode obtained by linear coupling of photonic and phononic modes. We described how hybrid photon-phonon blockade can be generated and detected in a driven nonlinear optomechanical superconducting system. Specifically, we considered the system composed of linearly coupled microwave and mechanical resonators with a superconducting qubit inserted in one of them.
We studied boson-number correlations in the photon, phonon, and hybrid modes in the system. By analysing steady-state second-order correlation functions, we found such parameter regimes of the system for which four different types of boson blockade and/or boson-induced tunnelling can be observed. Thus, we showed that bosons generated in the studied system can exhibit the sub-Poissonian (or super-Poissonian) boson-number statistics accompanied by boson antibunching in some cases or bunching in others. These results can be interpreted as four different types of blockade or tunnelling effects, as summarised in Table 2.
By tuning the pump frequency with respect to the energy levels of the hybrid system, which is driven via the SMR, we showed that it is possible to observe PB and PIT that can be explained by a large energy-level anharmonicity in the strong-coupling (or large-nonlinearity) regime. However, the time evolution of the second-order correlation function g (2) (τ ) oscillates due to the coupling g between the SMR and qubit as well as the hopping f between the SMR and QD. We showed that it is possible to induce PB in the hybrid mode c that survives for much longer delay times by driving the QD instead of the SMR.
We also predicted unconventional PB in the three modes in the weak-coupling (or weak-nonlinearity) regime using a non-Hermitian Hamiltonian approach based on neglecting quantum jumps. Our analytical approximate predictions are in a relatively good agreement with our precise master-equation solutions (including quantum jumps).
Moreover, as summarised in Table 1, we showed the possibility to observe eight different combinations of either PB or PIT in the three modes (a, b, and c) in different coupling regimes of this system. Thus, in particular, we found that the tunnelling effects in the photonic and phononic modes can lead, by their simple linear mixing, to the hybrid photon-phonon blockade effect.
Finally, we discussed two methods of detecting hybrid-mode correlations. One of them is based on measuring various moments of photons and phonons via balanced homodyne correlation measurements. While the other method is based on converting phonons of the hybrid mode into photons, by using a linear coupler acting as a multi-level SWAP gate, and then applying the standard optical HBT interferometry.
We believe that our study of the interplay between photons and phonons can lead to developing new experimental methods for controlling and testing the quantum states of mechanical systems with atom-cavitymechanics polaritons. We hope that our work can also stimulate research on quantum engineering with hybrid photon-phonon modes. , which corresponds to the super-Poissonian statistics of orders k = 2 , 3, and 4, which might be interpreted, as the induced tunnelling by one, two, and three photons. However, we can also find intermediate four out of six cases, which can be interpreted as non-standard types single-PB and/or single-PIT, and in some cases can be identified as multi-PB 11,26,29,30,53 . However, a detailed classification of such multi-PB and their interpretation is not at the focus of this paper. The presented results show only the possibility of generating in our system a plethora of various photon-phonon correlation effects, which can be revealed by higher-order correlation functions for the experimentally feasible parameters.

Analytical approach via non-Hermitian Hamiltonian in Eq. (21).
Here, we follow the method of Ref. 45 to derive the coefficients C n,m,k and C ′ n,m,k for n, m ∈ 0, 1, 2 and k = e, g , which appear in Eqs. (22) and (23). First we recall that the balanced linear coupler (or a balanced beam splitter) transformation, which leads to Eq. (24), if applied to the input Fock states |n a , n b � for n a + n b ≤ 2 yields: So, for the input state |� abq (t)� , given in Eq. (20), the output state of the balanced linear coupler can be represented as follows: where the superposition coefficients are: (31) g 234 = sgn log g (2) z (0), sgn log g (3) z (0), sgn log g (4) z (0) .
(33) |� cdq (t)� = C 00g |00g� + e −iω d t C 00e |00e� + C ′ 10g |10g� + C ′ 01g |01g� Table 3. Different predictions of the nth-order super-and sub-Poissonian statistics with n = 2, 3, 4 for the photon ( z = a ), phonon (b), and hybrid photon-phonon (c) modes, where g 234 is defined in Eq. (31). The cases marked with √ can be identified under both (1) nonresonance conditions, as shown in Fig. 14b, and (2) resonance conditions, as shown in Fig. 14a, except the case marked with * . www.nature.com/scientificreports/ We can calculate the coefficients C n a ,n b ,g/e iteratively 45 . For a single excitation and assuming the resonance case SMR = m = q = and κ a = κ b = κ , the steady-state superposition coefficients can be calculated from: where η = η b , � = ω i − ω p and ω SMR = ω m = ω q = ω . Moreover, we assume the weak-driving regime. So, in the first iteration, the contributions from the states with more than a single excitation, such as C 01e , C 11g , ..., are negligible. From Eq. (35), by comparing the coefficients with a single excitation, we can see that C 10g and C 00e are much larger than C 01g , because of a weak-pump amplitude η , and they can be written as In the second iteration, to include states with two excitations in total, the steady-state coefficients can be calculated from: where � κ = � − iκ/2 and � γ = � − iγ /2 . As can be seen from Eq. (37), we have So, to minimise C 02g , the minimalization of C 11g and C 01g is also required. Destructive interference between the direct and indirect excitation paths in the energy ladders of the total system can enable us minimising C 02g . This explains the occurrence of the dip in g (2) b (0) in the mode b, as a signature of PB. As clearly seen in Fig. 12a, the optimal PB in this mode occurs at � SMR /g = ±1.2 . The above equations lead us to analytical optimal conditions for the system parameters to maximise the sub-Poissonian character of the QD mode and, thus, to optimise the parameters for observing PB in the mode b. Given Eq. (39) for a single excitation and Eq. (41) for two excitations, which are calculated from Eq. (37), we show that the second-order correlation function calculated by this method and the master equation method both give very similar predictions, as shown in Fig. 12, where the thick curves are calculated based on the non-Hermitian Hamiltonian approach and the thin curves correspond to the master-equation approach for the modes a, b, and c.